function b_SAAD=Clogit_SAAD(Y_SAAD,X_SAAD,Z_SAAD)

choice=2:13;
ns=numel(choice);

Y=sortrows(Y_SAAD,[1,2]);
Y=Y(:,3:end);
Y=Y-1;  %!!!### UPDATED on 06/03/2020: Y MUST be from 1 to J!!!

X=sortrows(X_SAAD,[1,2]);
X=[ones(size(X,1),1),X(:,3:end)];

Z_SAAD=sortrows(Z_SAAD,[1,2,3]);
Z_SAAD(Z_SAAD(:,3)==1,:)=[];
Z=zeros(size(Y,1),size(Z_SAAD,2),ns);
for i=4:size(Z_SAAD,2)
    Z(:,i,:)=reshape(Z_SAAD(:,i),ns,[])';
end
Z(:,1:3,:)=[];


% X coefficients: constant bdtot profit teaching HHI totpop65  
b(:,3) = [-1.628; -0.000141; -0.0391; -0.118; 0.629; 0.238];
b(:,4) = [-0.182; -0.0135; 1.732; -14.05; 1.657; 0.238];
b(:,5) = [2.235; -0.0198; 0.258; -0.642; 0.484; 0.0498];
b(:,6) = [-0.392; 0.000626; -0.839; -0.109; 0.00824; 0.0704];
b(:,7) = [-2.492; 0.00123; -1.047; 0.217; 0.986; 0.238];
b(:,8) = [-13.85; 0.00218; -13.31; -0.761; 4.385; 0.985 ];
b(:,9) = [3.321; -0.0134; 1.645; -13.88; -0.185; -0.111];
b(:,10) = [-0.0236; 0.00123; -1.290; -1.024; 0.108; 0.0571];
b(:,11) = [-2.989; 0.00111; -0.580; -0.897; 1.164; 0.279];
b(:,12) = [-0.0234; -0.00184; -0.768; -1.265; 1.173; 0.190];
b(:,13) = [-5.196; -0.00152; 0.0493; -0.187; 1.926; 0.444];


% Z coefficients: vmsh*25%Q, vmsh*25-50%Q, vmsh50-75%Q, vmsh75%Q, vlag*25%Q, vlag*25-50%Q, vlag50-75%Q, vlag75%Q
bz=[-0.195; 0.830; 0.659; 0.433; 5.595; 5.858; 5.467; 5.089];

bOri=[reshape(b(:,3:end),[],1);bz];
%bAns=zeros(size(bOri));

options=optimset('MaxFunEvals',200000000000,'MaxIter',1500000000000,'TolX',1e-16,'Tolfun',1e-16,'GradObj','on','DerivativeCheck','off','FinDiffType','central');
startval = bOri;
baseAlt=1;
restrMat=[];

b_SAAD = fminunc('clogit',startval,options,restrMat,Y,X,Z,baseAlt);








